!================================================================================!
! This file is part of crest.
!
! Copyright (C) 2021 - 2023 Philipp Pracht
!
! crest is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! crest is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with crest.  If not, see <https://www.gnu.org/licenses/>.
!================================================================================!

module crest_calculator
!>--- types and readers
  use iso_fortran_env,only:wp => real64,int64
  use strucrd
  use calc_type
  use crest_calculator_printout
!>--- potentials and API's
  use subprocess_engrad !> driver exports for subprocesses
  use api_engrad  !> contains many potentials
!>--- other
  use constraints
  use nonadiabatic_module
  use lwoniom_module
!$ use omp_lib
  implicit none
!=========================================================================================!
!>--- private module variables and parameters
  private

!=========================================================================================!
!>--- RE-EXPORTS of module calc_type
  public :: calcdata              !> calculator main object
  public :: calculation_settings  !> different calculation objects (levels) within calcdata
  public :: jobtype               !> calculation type ID's
  public :: get_dipoles
!>--- RE-EXPORT of constraints
  public :: constraint
  public :: scantype
  public :: calc_constraint
!>--- RE-EXPORT of printout routines
  public :: calculation_summary
!=========================================================================================!

!>--- global engrad call counter
  integer(int64),public :: engrad_total = 0

!>--- public module routines
  public :: potential_core
  public :: engrad
  interface engrad
    module procedure :: engrad_mol
  end interface engrad
  public :: numgrad

  public :: calc_eprint
  interface calc_eprint
    module procedure :: calc_print_energies
  end interface calc_eprint

  public :: numhess1,numhess2
  public :: constrhess

!========================================================================================!
!========================================================================================!
contains  !> MODULE PROCEDURES START HERE
!========================================================================================!
!========================================================================================!

  subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
!***************************************************************
!* subroutine engrad
!* Main routine to perform some energy and gradient calculation
!*
!* If anything needs to be allocated we must avoid doing it
!* locally and do it in calcdata instead. This is to avoid
!* overhead generated by the allocation, which can be massive
!* for thousands of energy+gradient calls.
!***************************************************************
    implicit none
    type(coord),target :: mol
    type(calcdata) :: calc
    real(wp),intent(inout) :: energy
    real(wp),intent(inout) :: gradient(3,mol%nat)
    integer,intent(out) :: iostatus
    integer :: i,j,k,l,n,io,nocc
    real(wp) :: dum1,dum2
    real(wp) :: efix,gnorm
    type(coord),pointer :: molptr
    integer :: pnat
    logical :: useONIOM

!>--- reset
    call initsignal()
    energy = 0.0_wp
    gradient(:,:) = 0.0_wp

!**********************************
!>--- check for sane input
!**********************************
    dum1 = sum(mol%xyz)
    if (dum1 .ne. dum1) then !> NaN catch, we don't want to calculate garbage.
      iostatus = 1           !> For some builds I found this necessary because
      return                 !> OpenMP can get picky...
    end if

!>--- Calculation setup
    n = calc%ncalculations

    !$omp critical
    if (n > 0) then
      if (.not.allocated(calc%etmp)) allocate (calc%etmp(n),source=0.0_wp)
      if (.not.allocated(calc%grdtmp)) allocate (calc%grdtmp(3,mol%nat,n),source=0.0_wp)
      if (.not.allocated(calc%eweight)) then
        allocate (calc%eweight(n),source=0.0_wp)
        do i = 1,n
          calc%eweight(i) = calc%calcs(i)%weight
        end do
      end if
      !>--- count the engrad call
      engrad_total = engrad_total+1
    end if

!>--- update ONIOM geometries
    useONIOM = allocated(calc%ONIOM)
    if (useONIOM) then
      if (.not.allocated(calc%ONIOMmols)) then
        allocate (calc%ONIOMmols(calc%ONIOM%ncalcs))
      end if
      call ONIOM_update_geo(calc%ONIOM,mol,calc%ONIOMmols,calc%ONIOMmap)
    end if

    iostatus = 0
    dum1 = 1.0_wp
    dum2 = 1.0_wp
    if (n > 0) then
      calc%etmp = 0.0_wp
      calc%grdtmp = 0.0_wp
    end if
    !$omp end critical

!**************************************
!>--- Calculation of potentials
!**************************************
    if (n > 0) then
      if (calc%id > 0.and.calc%id > n) error stop 'Invalid calculator setup'

!==================================================================================!
!>--- loop over all calculations to be done
      do i = 1,calc%ncalculations

        !> skip through calculations we do not want
        if (calc%calcs(i)%refine_lvl /= calc%refine_stage) cycle

        !> Assign the molecule (necessary for ONIOM stuff)
        if (calc%calcs(i)%ONIOM_id /= 0) then
          !j = calc%calcs(i)%ONIOM_id
          j = calc%ONIOMrevmap(i)
          call ONIOM_associate_mol(calc%ONIOMmols(j),molptr)
        else
          molptr => mol
        end if
        pnat = molptr%nat

        !> also skip through if only one level was requested
        if (calc%id > 0.and.i .ne. calc%id.and..not.useONIOM) cycle
        if (.not.calc%calcs(i)%active) cycle

        !==========================================!
        !> the actual potential call
        !==========================================!
        call potential_core(molptr,calc,i,iostatus)
        !==========================================!
        !==========================================!

        if (iostatus /= 0) then
          return
        end if
      end do
!==================================================================================!

!***************************************************
!>--- Select energy and gradient construction
!***************************************************
      !>--- for ONIOM calculations, copy gradients to right positions
      !>--- and project with Jacobian
      !$omp critical
      if (useONIOM) then
        call calc_ONIOM_projection(calc)
      end if
      !$omp end critical

      select case (calc%id)
      case (0) !> the DEFAULT
        !>--- an option to add multiple energies and gradients accodring to weights
        !>--- which might be useful for additive contributions
        call calc_add_weighted_egrd(n,calc%eweight,calc%etmp,calc%grdtmp, &
        &                energy,gradient)

        !$omp critical
        if (useONIOM) then
          call ONIOM_engrad(calc%ONIOM,mol,energy,gradient)
        end if
        !$omp end critical

      case (1:)
        !>--- if calc%id is a positive integer, take e+grd from
        !>--- the respective specified calculation
        j = calc%id
        if (j <= calc%ncalculations) then
          if (useONIOM.and.calc%calcs(j)%ONIOM_id /= 0) then
            k = calc%calcs(j)%ONIOM_id
            l = calc%calcs(j)%ONIOM_highlowroot
            call ONIOM_get_fraggrad(calc%ONIOM,k,gradient,l,energy)
          else
            energy = calc%etmp(j)
            gradient = calc%grdtmp(:,:,j)
          end if
        end if

      case (-1)
        if (useONIOM) then
          write (stdout,'(a)') '**ERROR** MECP routines currently incompatible with ONIOM'
          error stop
        end if
        !>--- if calc%id is equal to -1, we are using the MECP routines:
        !>--- take non-adiabatic arithmetic mean (of first two calculations)
        call engrad_mean(mol%nat,calc%ncalculations,calc%etmp,calc%grdtmp, &
        &                energy,gradient)

      case default
        !>--- any other case does not return any energy and gradient
        !>--- from this part  of the calculator.
        !>--- constraints (see below) will be applied, however.
      end select

      gnorm = norm2(gradient)

      !>--- printout (to file or stdout)
      call calc_eprint(calc,energy,calc%etmp,gnorm)
    end if

!**********************************************
!>--- Constraints
!**********************************************
    if (calc%nconstraints > 0) then
      !$omp critical
      if (.not.allocated(calc%grdfix)) then
        allocate (calc%grdfix(3,mol%nat),source=0.0_wp)
      end if
      !$omp end critical

      !!$omp parallel &
      !!$omp shared(calc%cons,mol,n,energy,gradient) &
      !!$omp private(efix,i,calc%grdfix)
      do i = 1,calc%nconstraints
        efix = 0.0_wp
        calc%grdfix = 0.0_wp
        if (calc%nfreeze > 0) then
          call calc%cons(i)%addfreeze(calc%freezelist)
        end if
        if (calc%cons(i)%type >= 0) then
          !>--- structural constraints
          call calc_constraint(mol%nat,mol%xyz,calc%cons(i),efix,calc%grdfix)

        else if (allocated(calc%etmp).and.allocated(calc%grdtmp)) then
          !>--- non-adiabatic constraints
          if (n > 1) then
            call calc_nonadiabatic_constraint(mol%nat,calc%cons(i),n,calc%etmp,calc%grdtmp, &
            &    efix,calc%grdfix)
          else !> this "else" is necessary for constrained model hessians
            efix = 0.0_wp
            calc%grdfix = 0.0_wp
          end if
        end if
        !!$omp critical
        energy = energy+efix
        gradient = gradient+calc%grdfix
        !!$omp end critical
      end do
      !!$omp end parallel
    end if

!**********************************************
!>--- Frozen atoms
!**********************************************
    call calc%freezegrad(gradient)

    return
  end subroutine engrad_mol

!========================================================================================!
!========================================================================================!
!========================================================================================!

  subroutine potential_core(molptr,calc,id,iostatus)
!***************************************************************
!* subroutine porential_core
!*
!* wrapped main potential calls.
!* enrgies and gradients are saved to calc
!* id is used to select the calculation
!***************************************************************
    implicit none
    type(coord),intent(in)       :: molptr
    type(calcdata),intent(inout) :: calc
    integer,intent(in)           :: id
    integer,intent(out)          :: iostatus
    integer :: pnat
    real(wp) :: dum1,dum2

    if (id > calc%ncalculations) return

    pnat = molptr%nat

    !>--- select the calculation type
    select case (calc%calcs(id)%id)

    case (jobtype%xtbsys)
      !> xtb system call
      call xtb_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%generic)
      !> generic script/program call
      call generic_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%tblite)
      !> tblite api call
      call tblite_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%gfn0)
      !> GFN0-xTB api
      call gfn0_engrad(molptr,calc%calcs(id),calc%calcs(id)%g0calc,calc%etmp(id), &
      &                calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%gfn0occ)
      !> Special GFN0-xTB api given orbital population
      !> note the use of calc%g0calc instead of calc%calcs(id)%g0calc !
      call gfn0occ_engrad(molptr,calc%calcs(id),calc%g0calc,calc%etmp(id), &
      &                   calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%gfnff)
      !> GFN-FF api
      call gfnff_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%xhcff)
      !> XHCFF-lib
      call xhcff_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%turbomole)
      !> Turbomole-style SPs
      call turbom_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%orca)
      !> ORCA-style SPs
      call orca_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)

    case (jobtype%lj)
      !> Lennard-Jones potential calculation
      if (allocated(calc%calcs(id)%other)) then
        read (calc%calcs(id)%other,*) dum1,dum2
      else
        dum2 = 3.4_wp !> Argon-LJ distance (σ) in Angstroem
        dum1 = 0.238129_wp/627.50947428_wp !> Argon-LJ energy (ε) in atomic units
      end if
      call lj_engrad(molptr%nat,molptr%xyz*autoaa,dum1,dum2, &
      &              calc%etmp(id),calc%grdtmp(:,1:pnat,id))
      calc%grdtmp(:,:,id) = calc%grdtmp(:,:,id)*autoaa

    case default
      calc%etmp(id) = 0.0_wp
      calc%grdtmp(:,:,id) = 0.0_wp
    end select

  end subroutine potential_core

!========================================================================================!
!========================================================================================!
!========================================================================================!

  subroutine numgrad(mol,calc,angrad,numgradient)
!*******************************************************
!* subroutine numgrad
!* routine to perform a numerical gradient calculation
!*******************************************************
    implicit none
    type(coord),intent(inout) :: mol
    type(calcdata),intent(inout) :: calc
    real(wp),intent(in),optional :: angrad(3,mol%nat)
    real(wp),allocatable,intent(out),optional :: numgradient(:,:)
    integer :: i,j,k,l,ich,och,io

    real(wp) :: energy,el,er
    real(wp),allocatable :: grad(:,:)
    real(wp),allocatable :: ngrd(:,:)
    real(wp),parameter :: step = 0.0005_wp
    real(wp),parameter :: step2 = 0.5_wp/step

    allocate (grad(3,mol%nat),source=0.0_wp)
    allocate (ngrd(3,mol%nat),source=0.0_wp)

    do i = 1,mol%nat
      do j = 1,3
        call calc%dealloc_params()
        mol%xyz(j,i) = mol%xyz(j,i)+step
        call engrad(mol,calc,er,grad,io)

        call calc%dealloc_params()
        mol%xyz(j,i) = mol%xyz(j,i)-2*step
        call engrad(mol,calc,el,grad,io)

        mol%xyz(j,i) = mol%xyz(j,i)+step
        ngrd(j,i) = step2*(er-el)
      end do
    end do

    if (present(angrad)) then
      write (stdout,*)
      write (stdout,*) 'Numerical Gradient:'
      do i = 1,mol%nat
        write (*,'(3f18.8)') ngrd(1:3,i)
      end do

      write (stdout,*)
      write (stdout,*) 'Gradient Difference:'
      do i = 1,mol%nat
        write (stdout,'(3f18.8)') ngrd(1:3,i)-angrad(1:3,i)
      end do
    end if

    if (present(numgradient)) then
      call move_alloc(ngrd,numgradient)
    else
      deallocate (ngrd)
    end if
    deallocate (grad)

    return
  end subroutine numgrad

!========================================================================================!
  subroutine numhess1(nat,at,xyz,calc,hess,io)
!***************************************************
!* Calculate and return the Hessian for the
!* COMBINED energy and gradient from engrad,
!* including constraints
!***************************************************
    implicit none

    integer,intent(in) :: nat
    integer,intent(in) :: at(nat)
    real(wp),intent(in) :: xyz(3,nat)
    type(calcdata) :: calc
    real(wp),intent(out) :: hess(nat*3,nat*3)
    integer,intent(out)  :: io

    type(coord) :: mol !> coord type, so that the input remains unchanged
    real(wp) :: energy,el,er,hij
    real(wp),allocatable :: gradr(:,:),gradl(:,:)
    real(wp),allocatable :: gradr_tmp(:,:,:),gradl_tmp(:,:,:)
    real(wp),parameter :: step = 0.005_wp,step2 = 0.5_wp/step !0.00001_wp
    integer :: i,j,k,l,m,ii,jj

    hess = 0.0_wp
    io = 0
    mol%nat = nat
    mol%at = at
    mol%xyz = xyz

    allocate (gradr(3,mol%nat),source=0.0_wp) !dummy
    allocate (gradl(3,mol%nat),source=0.0_wp) !dummy

    do i = 1,mol%nat
      do j = 1,3
        ii = (i-1)*3+j
        gradr = 0.0_wp
        mol%xyz(j,i) = mol%xyz(j,i)+step
        call engrad(mol,calc,er,gradr,io)

        gradl = 0.0_wp
        mol%xyz(j,i) = mol%xyz(j,i)-2.0_wp*step
        call engrad(mol,calc,el,gradl,io)

        mol%xyz(j,i) = mol%xyz(j,i)+step

        do k = 1,mol%nat
          do l = 1,3
            jj = (k-1)*3+l
            hess(jj,ii) = (gradr(l,k)-gradl(l,k))*step2
          end do
        end do

      end do
    end do

    !Symmetrize Hessian
    do i = 1,nat*3
      do j = i,nat*3
        hij = (hess(i,j)+hess(j,i))*0.5_wp
        hess(i,j) = hij
        hess(j,i) = hij
      end do
    end do

    call engrad(mol,calc,el,gradl,io) !>- to get the gradient of the non-displaced structure

    deallocate (gradl,gradr)
    call mol%deallocate()
    return
  end subroutine numhess1

!========================================================================================!
  subroutine numhess2(nat,at,xyz,calc,hess,io)
!***************************************************
!* Calculate and return the Hessian for EACH of the
!* calculation levels, but contributions from
!* constraints are ignored.
!***************************************************
    implicit none

    integer,intent(in) :: nat
    integer,intent(in) :: at(nat)
    real(wp),intent(in) :: xyz(3,nat)
    type(calcdata) :: calc
    real(wp),intent(out) :: hess(nat*3,nat*3,calc%ncalculations)
    integer,intent(out)  :: io

    type(coord) :: mol !> coord type, so that the input remains unchanged
    real(wp) :: energy,el,er,hij
    real(wp),allocatable :: gradr(:,:),gradl(:,:)
    real(wp),allocatable :: gradr_tmp(:,:,:),gradl_tmp(:,:,:)
    real(wp),parameter :: step = 0.005_wp,step2 = 0.5_wp/step !0.00001_wp
    integer :: i,j,k,l,m,ii,jj
    logical :: rddipoles
    real(wp),allocatable :: dipr(:,:),dipl(:,:)

    hess = 0.0_wp
    io = 0
    mol%nat = nat
    mol%at = at
    mol%xyz = xyz

    allocate (gradr(3,mol%nat),source=0.0_wp) !dummy
    allocate (gradl(3,mol%nat),source=0.0_wp) !dummy

    allocate (gradr_tmp(3,mol%nat,calc%ncalculations),source=0.0_wp)
    allocate (gradl_tmp(3,mol%nat,calc%ncalculations),source=0.0_wp)

!>--- prepare potential dipole gradient readout
    rddipoles = (any(calc%calcs(:)%rddip))
    if (rddipoles) then
      do m = 1,calc%ncalculations
        if (calc%calcs(m)%rddip) then
          calc%calcs(m)%rddipgrad = .true.
          if (allocated(calc%calcs(m)%dipgrad)) deallocate (calc%calcs(m)%dipgrad)
          allocate (calc%calcs(m)%dipgrad(3,nat*3))
        end if
      end do
    end if

!>--- actual numerical Hessian loop
    do i = 1,mol%nat
      do j = 1,3
        ii = (i-1)*3+j
        mol%xyz(j,i) = mol%xyz(j,i)+step
        call engrad(mol,calc,er,gradr,io)

        gradr_tmp = calc%grdtmp
        if (rddipoles) call get_dipoles(calc,dipr)

        mol%xyz(j,i) = mol%xyz(j,i)-2.0_wp*step
        call engrad(mol,calc,el,gradl,io)

        gradl_tmp = calc%grdtmp
        if (rddipoles) call get_dipoles(calc,dipl)

        mol%xyz(j,i) = mol%xyz(j,i)+step

!>--- construct numerical Hessians from analytical gradients
        do m = 1,calc%ncalculations
          do k = 1,mol%nat
            do l = 1,3
              jj = (k-1)*3+l
              hess(jj,ii,m) = (gradr_tmp(l,k,m)-gradl_tmp(l,k,m))*step2
            end do
          end do
        end do

!>--- create dipole gradients
        if (rddipoles) then
          do m = 1,calc%ncalculations
            if (calc%calcs(m)%rddip) then
              calc%calcs(m)%dipgrad(1:3,ii) = (dipr(1:3,m)-dipl(1:3,m))*step2
            end if
          end do
        end if

      end do
    end do

!>--- Symmetrize Hessian
    do m = 1,calc%ncalculations
      do i = 1,nat*3
        do j = i,nat*3
          hij = (hess(i,j,m)+hess(j,i,m))*0.5_wp
          hess(i,j,m) = hij
          hess(j,i,m) = hij
        end do
      end do
    end do

    call engrad(mol,calc,el,gradl,io) !>- to get the gradient of the non-displaced s

    if (allocated(dipl)) deallocate (dipl)
    if (allocated(dipr)) deallocate (dipr)
    deallocate (gradl_tmp,gradr_tmp)
    deallocate (gradl,gradr)
    call mol%deallocate()
    return
  end subroutine numhess2

!========================================================================================!
  subroutine constrhess(nat,at,xyz,calc,phess)
!*********************************************************
!* subroutine constrhess
!* routine to perform a numerical Hessian calculation
!* but ONLY include contributions of the constraints.
!*
!* phess is the packed Hessian on which the constraint
!* contributions are added.
!*********************************************************
    implicit none

    integer,intent(in) :: nat
    integer,intent(in) :: at(nat)
    real(wp),intent(in) :: xyz(3,nat)
    type(calcdata),intent(in) :: calc
    real(wp),intent(inout) :: phess((nat*3)*((nat*3)+1)/2)

    type(calcdata) :: dummycalc
    integer :: tc_backup
    real(wp) :: energy,el,er
    real(wp),allocatable :: hess(:,:)
    logical :: consgeo
    integer :: i,j,k,n3,io,ncons

    ncons = calc%nconstraints
    if (ncons <= 0) return
    !>--- skip if only nonadiabatic constraints

    !> new dummy calculation
    !> set to something arbitrary so that only constraints are considered
    !> and copy the neccesities
    dummycalc%id = -1000
    dummycalc%ncalculations = 0
    dummycalc%pr_energies = .false.

    !> copy the constraints
    dummycalc%nfreeze = calc%nfreeze
    dummycalc%nconstraints = ncons
    !$omp critical
    allocate (dummycalc%cons(ncons))
    do i = 1,ncons
      dummycalc%cons(i) = calc%cons(i)
    end do
    if (calc%nfreeze > 0) then
      dummycalc%freezelist = calc%freezelist
    end if
    n3 = nat*3
    allocate (hess(n3,n3),source=0.0_wp)
    !$omp end critical

    call numhess1(nat,at,xyz,dummycalc,hess,io)

    k = 0
    do i = 1,n3
      do j = 1,i
        k = k+1
        phess(k) = phess(k)+0.5_wp*(hess(j,i)+hess(i,j))
      end do
    end do

    deallocate (hess)
    return
  end subroutine constrhess

!==========================================================================================!
  subroutine calc_print_energies(calc,energy,energies,gnorm,chnl)
    implicit none
    type(calcdata),intent(inout) :: calc
    real(wp),intent(in) :: energy
    real(wp),intent(in) :: energies(:)
    real(wp),intent(in) :: gnorm
    integer,intent(in),optional :: chnl
    integer :: i,j,l
    character(len=20) :: atmp
    character(len=:),allocatable :: btmp
    logical :: ex

    if (.not.present(chnl)) then
      if (.not.calc%pr_energies) return
    end if
    btmp = ''
    write (atmp,'(f20.12)') energy
    btmp = btmp//atmp
    l = size(energies,1)
    do i = 1,l
      write (atmp,'(f20.12)') energies(i)
      btmp = btmp//atmp
    end do
    write (atmp,'(e20.5)') gnorm
    btmp = btmp//atmp
    if (present(chnl)) then
      write (chnl,'(a)') btmp
    else
      write (calc%eout_unit,'(a)') btmp
    end if
    deallocate (btmp)
    return
  end subroutine calc_print_energies

!==========================================================================================!
  subroutine calc_add_weighted_egrd(ncalc,weights,energies,gradients,e,grd)
!***************************************************************************
!* A subroutine that adds (further) energies and gradients to
!* e and grd, according to given weights.
!***************************************************************************
    implicit none
    !> INPUT
    integer,intent(in) :: ncalc
    real(wp),intent(in) :: weights(ncalc)
    real(wp),intent(in) :: energies(ncalc)
    real(wp),intent(in) :: gradients(:,:,:)
    !> OUTPUT
    real(wp),intent(inout) :: e
    real(wp),intent(inout) :: grd(:,:)
    !> LOCAL
    integer :: i
    do i = 1,ncalc
      e = e+weights(i)*energies(i)
      grd(:,:) = grd(:,:)+weights(i)*gradients(:,:,i)
    end do
  end subroutine calc_add_weighted_egrd

!==========================================================================================!

  subroutine calc_ONIOM_projection(calc)
!*******************************************
!* Iterate through the ONIOM data and
!* place the correct energies and gradients
!*******************************************
    implicit none
    type(calcdata),intent(inout) :: calc
    integer :: n,i,j,k,l,l1,l2
    integer :: natp,trunat
#ifdef WITH_LWONIOM
    if (allocated(calc%ONIOM).and.calc%ncalculations > 0) then
      trunat = maxval(calc%ONIOMmols(:)%nat)
      do i = 1,calc%ONIOM%nfrag
        l1 = calc%ONIOM%calcids(1,i)
        l2 = calc%ONIOM%calcids(2,i)
        natp = calc%ONIOMmols(l1)%nat

        call calc%ONIOM%fragment(i)%gradient_distribution(  &
        &    calc%etmp(l1),calc%grdtmp(1:3,1:natp,l1), &
        &    calc%etmp(l2),calc%grdtmp(1:3,1:natp,l2))
        call calc%ONIOM%fragment(i)%jacobian(trunat)
      end do
    end if
#endif
  end subroutine calc_ONIOM_projection

!==========================================================================================!
!==========================================================================================!
end module crest_calculator
